# This tutorial will use Monte Carlo methods to simulate the output of a
# Photomultiplier Tube (PMT) once triggered by the detection of one or more
# photons.
# This example follows notes posted online which you can find at the following
# url: http://superk.physics.sunysb.edu/~mcgrew/phy310/lectures/phy310-lecture-06-2007.pdf
import numpy as np
import matplotlib.pyplot as plt
import statistics
# In all Monte Carlo simulations it is necessary to generate random or
# pseudo-random numbers. The Python command 'np.random.uniform()'
# from the NumPy module generates random numbers uniformly distributed
# between zero and one.
np.random.uniform()
0.7961852081589491
# First, suppose N photons hit the photocathode. Determine the number of
# photoelectrons (pe) that are generated. Assume that the incoming light is
# 400 nm and that the quantum efficiency of the photocathode is 0.23.
#(Following http://superk.physics.sunysb.edu/~mcgrew/phy310/lectures/phy310-lecture-06-2007.pdf)
N = 2363
QE = 0.23
pe = 0
for i in range(N):
if np.random.uniform() < QE:
pe += 1
pe
546
# When an electron with energy E hits a Dynode the average number of
# secondary electrons liberated is alpha*E (i.e. number of secondary electrons
# is proportional to the energy of the incoming electron) where alpha is a
# contant. Since we are counting electrons, the distribution of liberated
# electrons follows the Poisson distribution. In Python random numbers
# drawn from a Poisson distribution with mean mu are generated using
# poissrnd(mu). As an example, below we generate random integers drawn from
# a Poisson parent distribtuion with mean 2000.
mu = 2000
Ydist = []
for i in range(1000):
Ydist = Ydist + [np.random.poisson(mu)]
print('mean:', statistics.mean(Ydist))
print('variance:', statistics.stdev(Ydist)**2)
nbins = 25
plt.hist(Ydist, nbins, color = 'c', edgecolor = 'k')
plt.xlabel('Y')
plt.ylabel('counts');
mean: 1998.539 variance: 1942.5470260260263
# A photomultiplier tube (PMT) consists of a photocathode followed by a
# series of dynodes maintained at different electric potentials and then
# finally an anode. When a single photon is incident on the photocathode
# it either produces a photoelectron or it doesn't. The probability that
# it produces a photoelectron is determined by the quantum efficiency
# QE of the photocathode. For this exercise we will assume that QE = 0.23.
# If a photoelectron is produced, it is accelerated towards the first dynode by
# means of a potential difference. We assume that all electrons, whether
# produced at the photocathode or one of the dynodes, start with zero kinetic
# energy. Therefore, the energy an electron gains is simply its charge times
# the potential difference between its starting and final positions. When
# an electron collides with a dynode, secondary electrons are produced. The
# average number of secondary electrons produced is proportional to the energy
# of the incoming electron. In this problem, we assume that an electron
# accelerated through 20 V will, on average, produce one secondary electron
# upon colliding with the dynode. Because we are "counting" electrons, the
# distribution of secondary electrons generated will follow a Poisson
# distribution.
QE = 0.23; # Set the quantum efficeincy of the photocathode
maxi = int(1e4) # Set the number of Monte Carlo iterations
dynode = [0, 150, 300, 450, 600, 750, 850] # Set the number of dynodes and
# their voltages
dist = []
from datetime import datetime
print(datetime.now()) # Print the time at the simulation was started
for i in range(maxi): # This loop runs the Monte Carlo simulation maxi times
# (10,000 in this case)
if np.random.uniform() < QE: # Generate a random number between [0, 1] and check to see
# if it's greater than the quantum efficiency
electrons = 1 # If true, then a photoelectron is produced
for j in range(len(dynode)-1): # This loop will examine the secondary
# electrons produced at each of the dynodes
mu = electrons*(dynode[j+1]-dynode[j])/20 # Mean number of
# electrons generated at a dynode determined from the number of
# incoming electons and the energy of the incoming electrons.
# One electron accelerated through 20 V will on anverage produce
# one electron
electrons = electrons + np.random.poisson(mu)
dist = dist + [electrons] # Add the total number of electrons detected at
# the anode for this trial to a 'dist' list.
print(datetime.now())
2021-03-11 11:22:05.176837 2021-03-11 11:22:05.230722
print('The fraction of trials that produced photoelectrons is', len(dist)/maxi)
# Determines what fraction of maxi trials actually
# produced electrons at the PMT anode. It should come out to be close the
# QE = 0.23.
The fraction of trials that produced photoelectrons is 0.2286
import statistics
print('The mean of the distribution was', statistics.mean(dist))
# Calculate the mean of the distribution
xbins = np.linspace(0, 2e6, 101) # Set bin ranges for the histogram.
# 100 equally-sized bins fron 0 to 2e6.
plt.hist(dist, xbins, color = 'c', edgecolor = 'k') # Plot the histgram of the
# simulated number of anode electrons
plt.xlim((0, 2e6))
plt.ylim((0, 500))
plt.xlabel('Number of Anode Electons')
plt.ylabel('Counts');
The mean of the distribution was 268264.28783902014
counts, edges = np.histogram(dist, xbins) # Here we use 'np.histogram' to have
# Python tell us how many counts there are in each of our 100 bins.
counts
array([ 0, 1, 5, 16, 21, 43, 76, 97, 146, 140, 185, 198, 202, 200, 169, 166, 147, 103, 92, 80, 51, 48, 30, 20, 22, 11, 6, 2, 4, 3, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int64)
# Calculate the centre positions of the bins from the edge locations.
spacing = (edges[1] - edges[0])/2
centres = edges[1:] - spacing
centres
array([ 10000., 30000., 50000., 70000., 90000., 110000., 130000., 150000., 170000., 190000., 210000., 230000., 250000., 270000., 290000., 310000., 330000., 350000., 370000., 390000., 410000., 430000., 450000., 470000., 490000., 510000., 530000., 550000., 570000., 590000., 610000., 630000., 650000., 670000., 690000., 710000., 730000., 750000., 770000., 790000., 810000., 830000., 850000., 870000., 890000., 910000., 930000., 950000., 970000., 990000., 1010000., 1030000., 1050000., 1070000., 1090000., 1110000., 1130000., 1150000., 1170000., 1190000., 1210000., 1230000., 1250000., 1270000., 1290000., 1310000., 1330000., 1350000., 1370000., 1390000., 1410000., 1430000., 1450000., 1470000., 1490000., 1510000., 1530000., 1550000., 1570000., 1590000., 1610000., 1630000., 1650000., 1670000., 1690000., 1710000., 1730000., 1750000., 1770000., 1790000., 1810000., 1830000., 1850000., 1870000., 1890000., 1910000., 1930000., 1950000., 1970000., 1990000.])
plt.semilogy(centres, counts, 'bo', fillstyle = 'none') # Make a semi-log plot
# of the counts versus the number of detected anode electrons
plt.xlim((0, 2e6))
plt.ylim((.8, 500))
plt.xlabel('Number of Anode Electrons')
plt.ylabel('Counts')
Text(0, 0.5, 'Counts')
# This next block of code is very similar to the code presented above. This time,
# however, we imagine that there are 8 photons incident on the photocathode
# instead of one. Now there can be 0, 1, 2, ..., or 8 photoelectrons
# generated and then accelerated towards the first dynode. How does the
# distribution of electrons arriving at the anode change? We only require
# a relatively minor change to the code to study this problem. Plot the
# histograms using the same scale and the same binwidths to make comparisons
# with the previous results easy.
QE = 0.23 # Set the quantum efficeincy of the photocathode
maxi = int(1e4) # Set the number of Monte Carlo iterations
dynode = [0, 150, 300, 450, 600, 750, 850] # Set the number of dynodes and
# their voltages
numPhotons = 8 # Set the number of photons incident on the photocathode.
dist = []
print(datetime.now()) # Print the time at the simulation was started
for i in range(maxi): # This loop runs the Monte Carlo simulation maxi times
# (10,000 in this case)
pe = 0 # Set the initial number of photoelectrons to zero.
for k in range(numPhotons): # This loop individually steps through each of the
# photons incident on the photocathode
if np.random.uniform() < QE: # Generate a random number between [0, 1] and check to see
# if it's greater than the quantum efficiency
pe += 1 # Increment pe by one every time a photoelectron is
# produced
if pe > 0: # If there are any photoelectrons, then do something. If there
# are not photoelectrons then do nothing
electrons = pe # Set the initial number of electrons equal to the
# number of photoelectrons
for j in range(len(dynode)-1): # This loop will examine the secondary
# electrons produced at each of the dynodes
mu = electrons*(dynode[j+1]-dynode[j])/20; # Mean number of
# electrons generated at a dynode determined from the number of
# incoming electons and the energy of the incoming electrons.
# One electron accelerated through 20 V will on anverage produce
# one electron
electrons = electrons + np.random.poisson(mu)
dist = dist + [electrons] # Add the total number of electrons detected at
# the anode for this trial to a 'dist' list.
print(datetime.now()) # Print the time at the simulation completed.
print('The fraction of trials that produced photoelectrons is', len(dist)/maxi)
# Determines what fraction of maxi trials actually
# produced electrons at the PMT anode.
print('The mean of the distribution was', statistics.mean(dist))
# Calculate the mean of the distribution
plt.figure() # Generate a new figure
xbins = np.linspace(0, 2e6, 101) # Set bin ranges for the histogram.
# 100 equally-sized bins fron 0 to 2e6.
plt.hist(dist, xbins, color = 'c', edgecolor = 'k') # Plot the histgram of the
# simulated number of anode electrons
plt.xlim((0, 2e6))
plt.ylim((0, 500))
plt.xlabel('Number of Anode Electons')
plt.ylabel('Counts')
counts, edges = np.histogram(dist, xbins) # Here we use 'np.histogram' to have
# Python tell us how many counts there are in each of our 100 bins.
# Calculate the centre positions of the bins from the edge locations.
spacing = (edges[1] - edges[0])/2
centres = edges[1:] - spacing
plt.figure() # Start a new figure.
plt.semilogy(centres, counts, 'bo', fillstyle = 'none') # Make a semi-log plot
# of the counts versus the number of detected anode electrons
plt.xlim((0, 2e6))
plt.ylim((.8, 500))
plt.xlabel('Number of Anode Electrons')
plt.ylabel('Counts');
2021-03-11 11:24:58.852798 2021-03-11 11:24:59.142357 The fraction of trials that produced photoelectrons is 0.8761 The mean of the distribution was 560754.4641022715
# Finally, this last block of code runs the same PMT simulation above four
# times, each time using a different number of incoming photons [1, 2, 4, 8].
# This chunk of code doesn't rely on anything that came before it. These 44 lines
# of code will simulate the expected distributions of anode electrons for 4
# difference sets of incoming photons. There real virtue of the Monte Carlo
# simulation is that we can now vary properties of the PMT with trivial
# modifications to the code below and systematically study the effects.
# For example, what if there were mode dynodes? All we have to do is modify
# the line dynode:=[0,150,300,450,600,750,800]. Alternatively, we could
# keep the number of dynodes fixed and modify the potential applied to the
# dynodes. Of course, after making this relatively simple simulation work,
# we could make modifications to make it more sophisticated. What if 8
# photons are directed towards the photocathode, but they arrive at slightly
# different times. What does the current pulse at the anode look like?
# That's not a problem that we'll tackle here, but it does demonstrate the
# versitility of the Monte Carlo method.
QE = 0.23 # Set the quantum efficeincy of the photocathode
maxi = int(1e4) # Set the number of Monte Carlo iterations
dynode = [0, 150, 300, 450, 600, 750, 850] # Set the number of dynodes and
# their voltages
numPhotons = [1, 2, 4, 8] # Set the number of photons incident on the
# photocathode. We will repeat the simulation for 4 different numners of
# incoming photons
print(datetime.now()) # Print the time at the simulation was started
eventsList = []
distList = []
for m in range (len(numPhotons)):
dist = []
plt.figure() # Generate a new figure
for i in range(maxi): # This loop runs the Monte Carlo simulation maxi times
# (10,000 in this case)
pe = 0 # Set the initial number of photoelectrons to zero.
for k in range (numPhotons[m]): # This loop individually steps through each of the
# photons incident on the photocathode
if np.random.uniform() < QE: # Generate a random number between [0, 1]
# and check to see if it's greater than the quantum efficiency
pe += 1 # Increment pe by one every time a photoelectron is
# produced
if pe > 0: # If there are any photoelectrons, then do something. If there
# are not photoelectrons then do nothing
electrons = pe # Set the initial number of electrons equal to the
# number of photoelectrons
for j in range(len(dynode)-1): # This loop will examine the secondary
# electrons produced at each of the dynodes
mu = electrons*(dynode[j+1]-dynode[j])/20 # Mean number of
# electrons generated at a dynode determined from the number of
# incoming electons and the energy of the incoming electrons.
# One electron accelerated through 20 V will on anverage produce
# one electron
electrons = electrons + np.random.poisson(mu)
dist = dist + [electrons] # Add the total number of electrons detected at
# the anode for this trial to a 'dist' list.
# This code can take some time to execute, so provide some intermediate
# outputs to track progress
print('Iteration', m)
print('The fraction of trials that produced photoelectrons is', len(dist)/maxi)
# Determines what fraction of maxi trials actually
# produced electrons at the PMT anode.
print('The mean of the distribution was', statistics.mean(dist))
# Calculate the mean of the distribution
print(datetime.now()) # Print the time that this iteration of the simulation completed
eventsList = eventsList + [[numPhotons[m], len(dist)/maxi]] # Keep track of
# how many sets of incident photons actually resulted in a signal at
# the anode.
distList = distList + [dist] # Store the distributions recorded into a set
xbins = xbins = np.linspace(0, 2e6, 101) # Set bin ranges for the histogram.
# 100 equally-sized bins fron 0 to 2e6.
plt.hist(dist, xbins, color = 'c', edgecolor = 'k')
#Plot the histgram of the simulated number of anode electrons
plt.xlim((0, 2e6))
plt.ylim((0, 500))
plt.xlabel('Number of Anode Electons')
plt.ylabel('Counts')
counts, edges = np.histogram(dist, xbins) # Here we use 'np.histogram' to have
# Python tell us how many counts there are in each of our 100 bins.
# Calculate the centre positions of the bins from the edge locations.
spacing = (edges[1] - edges[0])/2
centres = edges[1:] - spacing
plt.figure() # Start a new figure.
plt.semilogy(centres, counts, 'bo', fillstyle = 'none') # Make a semi-log plot
# of the counts versus the number of detected anode electrons
plt.xlim((0, 2e6))
plt.ylim((0.8, 500))
plt.xlabel('Number of Anode Electrons')
plt.ylabel('Counts');
2021-03-11 11:25:35.795103 Iteration 0 The fraction of trials that produced photoelectrons is 0.2217 The mean of the distribution was 262444.55435272894 2021-03-11 11:25:35.842940 Iteration 1 The fraction of trials that produced photoelectrons is 0.4072 The mean of the distribution was 302301.53806483303 2021-03-11 11:25:36.061257 Iteration 2 The fraction of trials that produced photoelectrons is 0.6407 The mean of the distribution was 378459.22553457157 2021-03-11 11:25:36.291688 Iteration 3 The fraction of trials that produced photoelectrons is 0.8822 The mean of the distribution was 560186.1229879846 2021-03-11 11:25:36.651354
# Using the data stored in 'eventsSet', generate a plot of the average number
# of electrons detected at the anode versus the number of photons incident
# on the photocathode.
y = []
for i in range(len(eventsList)):
y = y + [eventsList[i][1]]
plt.plot(numPhotons, y, 'ro', linewidth = 2, fillstyle = 'none')
plt.xlabel('Number of Incident Photons')
plt.ylabel('Average Number of Electrons at Anode');
# Finally, using the distributions stored in the set 'distSet', plot the
# simulated distributions for the different numbers of incident photons all on
# the same semi-log plot.
plt.figure(figsize=(12,10))
colours = ['ko', 'bo', 'ro', 'go']
for i in range(len(numPhotons)):
counts, edges = np.histogram(distList[i], xbins)
spacing = (edges[1] - edges[0])/2
centres = edges[1:] - spacing
plt.semilogy(centres, counts, colours[i], fillstyle = 'none', linewidth = 2)
plt.xlim((0, 2e6))
plt.ylim((0.8, 500))
plt.xlabel('Number of Anode Electrons')
plt.ylabel('Counts')
plt.legend(('1 photon', '2 photons', '4 photons', '8 photons'));